{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "7f7d69f9",
      "metadata": {
        "id": "7f7d69f9"
      },
      "source": [
        "Copyright © 2024 Gurobi Optimization, LLC\n",
        "\n",
        "# Price Optimization with Competing Products\n",
        "Finding the delicate balance between price and demand is a difficult problem to solve and one that is prevalent in a number of huge industries like retail, e-commerce, ticketing, and hospitality.\n",
        "\n",
        "It just so happens that this is a problem that data scientists have recently been getting better and better at addressing. Still, a key piece is missing: How can I make the right pricing decision given all the other constraints and business rules that exist for this problem? That's where optimization comes in.\n",
        "\n",
        "In this scenario, we have a product that has several categories, but we have a limited amount of \"space\" for our products. This could be shelf or warehouse space in retail, seats for events, or for airline ticketing, or rooms at a hotel.\n",
        "\n",
        "This problem considers two similar products that are offered where it's our job to use the data available with some optimization know-how to determine the optimal mix of products to offer that maximize revenue while also adhering to a few other business rules. First, we'll create a predictive model to forecast sales based on the prices of each product. Then we'll build an optimization model to find this optimal mix. Finally, we’ll also use the Gurobi-sponsored open-source package Gurobi Machine Learning to seamlessly combine the features of a machine learning model with a decision of an optimization model.\n",
        "\n",
        "Let’s get started!"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a5846a14",
      "metadata": {
        "id": "a5846a14"
      },
      "source": [
        "## Load required packages\n",
        "If you have a Gurobi license you can skip the installation of `gurobipy`, but always make sure you have the [latest version](https://www.gurobi.com/downloads/gurobi-software) available."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "1c3276f0",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1c3276f0",
        "outputId": "8b9efd0a-df20-42e8-99b7-28eb2ae2b3e5"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: gurobipy in /usr/local/lib/python3.10/dist-packages (11.0.0)\n"
          ]
        }
      ],
      "source": [
        "%pip install gurobipy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "3024295a",
      "metadata": {
        "id": "3024295a"
      },
      "outputs": [],
      "source": [
        "import gurobipy as gp\n",
        "from gurobipy import GRB\n",
        "\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import seaborn as sns\n",
        "import matplotlib.pyplot as plt\n",
        "import warnings\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn import tree\n",
        "\n",
        "warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d8c28cfd",
      "metadata": {
        "id": "d8c28cfd"
      },
      "source": [
        "## Start with some data analysis\n",
        "\n",
        "This data contains prices and sales for two of our competing products and was generated using another script, which can be found [here](). Let's load the data and take a quick look."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "eea57a4c",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "eea57a4c",
        "outputId": "28a70699-d016-4f52-96e3-bc9badf6e0b0"
      },
      "outputs": [
        {
          "data": {
            "application/vnd.google.colaboratory.intrinsic+json": {
              "summary": "{\n  \"name\": \"df\",\n  \"rows\": 1000,\n  \"fields\": [\n    {\n      \"column\": \"p[1]\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 14.409836071934338,\n        \"min\": 300.0,\n        \"max\": 400.0,\n        \"num_unique_values\": 911,\n        \"samples\": [\n          350.69,\n          359.04,\n          351.26\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"p[2]\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 29.093624863457553,\n        \"min\": 100.0,\n        \"max\": 300.0,\n        \"num_unique_values\": 947,\n        \"samples\": [\n          204.43,\n          211.29,\n          175.08\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"n[1]\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 29.65004858986544,\n        \"min\": 0.0,\n        \"max\": 200.0,\n        \"num_unique_values\": 148,\n        \"samples\": [\n          70.0,\n          105.0,\n          48.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    }\n  ]\n}",
              "type": "dataframe",
              "variable_name": "df"
            },
            "text/html": [
              "\n",
              "  <div id=\"df-78c5f739-5ddb-4a26-9e2e-3e411556b06c\" class=\"colab-df-container\">\n",
              "    <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>p[1]</th>\n",
              "      <th>p[2]</th>\n",
              "      <th>n[1]</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>356.12</td>\n",
              "      <td>197.67</td>\n",
              "      <td>108.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>358.05</td>\n",
              "      <td>189.68</td>\n",
              "      <td>66.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>340.79</td>\n",
              "      <td>260.35</td>\n",
              "      <td>130.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>353.76</td>\n",
              "      <td>133.53</td>\n",
              "      <td>55.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>341.37</td>\n",
              "      <td>229.80</td>\n",
              "      <td>91.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>...</th>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>995</th>\n",
              "      <td>357.63</td>\n",
              "      <td>241.54</td>\n",
              "      <td>68.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>996</th>\n",
              "      <td>352.58</td>\n",
              "      <td>212.95</td>\n",
              "      <td>87.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>997</th>\n",
              "      <td>355.28</td>\n",
              "      <td>189.50</td>\n",
              "      <td>94.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>998</th>\n",
              "      <td>369.75</td>\n",
              "      <td>166.33</td>\n",
              "      <td>51.0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>999</th>\n",
              "      <td>349.31</td>\n",
              "      <td>222.07</td>\n",
              "      <td>114.0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "<p>1000 rows × 3 columns</p>\n",
              "</div>\n",
              "    <div class=\"colab-df-buttons\">\n",
              "\n",
              "  <div class=\"colab-df-container\">\n",
              "    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-78c5f739-5ddb-4a26-9e2e-3e411556b06c')\"\n",
              "            title=\"Convert this dataframe to an interactive table.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n",
              "    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "\n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    .colab-df-buttons div {\n",
              "      margin-bottom: 4px;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "    <script>\n",
              "      const buttonEl =\n",
              "        document.querySelector('#df-78c5f739-5ddb-4a26-9e2e-3e411556b06c button.colab-df-convert');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      async function convertToInteractive(key) {\n",
              "        const element = document.querySelector('#df-78c5f739-5ddb-4a26-9e2e-3e411556b06c');\n",
              "        const dataTable =\n",
              "          await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                    [key], {});\n",
              "        if (!dataTable) return;\n",
              "\n",
              "        const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "          + ' to learn more about interactive tables.';\n",
              "        element.innerHTML = '';\n",
              "        dataTable['output_type'] = 'display_data';\n",
              "        await google.colab.output.renderOutput(dataTable, element);\n",
              "        const docLink = document.createElement('div');\n",
              "        docLink.innerHTML = docLinkHtml;\n",
              "        element.appendChild(docLink);\n",
              "      }\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "\n",
              "<div id=\"df-688f406a-bf03-42ab-8f80-78cfe2da5da0\">\n",
              "  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-688f406a-bf03-42ab-8f80-78cfe2da5da0')\"\n",
              "            title=\"Suggest charts\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "<svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "     width=\"24px\">\n",
              "    <g>\n",
              "        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n",
              "    </g>\n",
              "</svg>\n",
              "  </button>\n",
              "\n",
              "<style>\n",
              "  .colab-df-quickchart {\n",
              "      --bg-color: #E8F0FE;\n",
              "      --fill-color: #1967D2;\n",
              "      --hover-bg-color: #E2EBFA;\n",
              "      --hover-fill-color: #174EA6;\n",
              "      --disabled-fill-color: #AAA;\n",
              "      --disabled-bg-color: #DDD;\n",
              "  }\n",
              "\n",
              "  [theme=dark] .colab-df-quickchart {\n",
              "      --bg-color: #3B4455;\n",
              "      --fill-color: #D2E3FC;\n",
              "      --hover-bg-color: #434B5C;\n",
              "      --hover-fill-color: #FFFFFF;\n",
              "      --disabled-bg-color: #3B4455;\n",
              "      --disabled-fill-color: #666;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart {\n",
              "    background-color: var(--bg-color);\n",
              "    border: none;\n",
              "    border-radius: 50%;\n",
              "    cursor: pointer;\n",
              "    display: none;\n",
              "    fill: var(--fill-color);\n",
              "    height: 32px;\n",
              "    padding: 0;\n",
              "    width: 32px;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart:hover {\n",
              "    background-color: var(--hover-bg-color);\n",
              "    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "    fill: var(--button-hover-fill-color);\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart-complete:disabled,\n",
              "  .colab-df-quickchart-complete:disabled:hover {\n",
              "    background-color: var(--disabled-bg-color);\n",
              "    fill: var(--disabled-fill-color);\n",
              "    box-shadow: none;\n",
              "  }\n",
              "\n",
              "  .colab-df-spinner {\n",
              "    border: 2px solid var(--fill-color);\n",
              "    border-color: transparent;\n",
              "    border-bottom-color: var(--fill-color);\n",
              "    animation:\n",
              "      spin 1s steps(1) infinite;\n",
              "  }\n",
              "\n",
              "  @keyframes spin {\n",
              "    0% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "      border-left-color: var(--fill-color);\n",
              "    }\n",
              "    20% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    30% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    40% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    60% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    80% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "    90% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "  }\n",
              "</style>\n",
              "\n",
              "  <script>\n",
              "    async function quickchart(key) {\n",
              "      const quickchartButtonEl =\n",
              "        document.querySelector('#' + key + ' button');\n",
              "      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n",
              "      quickchartButtonEl.classList.add('colab-df-spinner');\n",
              "      try {\n",
              "        const charts = await google.colab.kernel.invokeFunction(\n",
              "            'suggestCharts', [key], {});\n",
              "      } catch (error) {\n",
              "        console.error('Error during call to suggestCharts:', error);\n",
              "      }\n",
              "      quickchartButtonEl.classList.remove('colab-df-spinner');\n",
              "      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n",
              "    }\n",
              "    (() => {\n",
              "      let quickchartButtonEl =\n",
              "        document.querySelector('#df-688f406a-bf03-42ab-8f80-78cfe2da5da0 button');\n",
              "      quickchartButtonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "    })();\n",
              "  </script>\n",
              "</div>\n",
              "\n",
              "  <div id=\"id_f05dd4cf-2dc8-426b-8b2b-c1d115c5d749\">\n",
              "    <style>\n",
              "      .colab-df-generate {\n",
              "        background-color: #E8F0FE;\n",
              "        border: none;\n",
              "        border-radius: 50%;\n",
              "        cursor: pointer;\n",
              "        display: none;\n",
              "        fill: #1967D2;\n",
              "        height: 32px;\n",
              "        padding: 0 0 0 0;\n",
              "        width: 32px;\n",
              "      }\n",
              "\n",
              "      .colab-df-generate:hover {\n",
              "        background-color: #E2EBFA;\n",
              "        box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "        fill: #174EA6;\n",
              "      }\n",
              "\n",
              "      [theme=dark] .colab-df-generate {\n",
              "        background-color: #3B4455;\n",
              "        fill: #D2E3FC;\n",
              "      }\n",
              "\n",
              "      [theme=dark] .colab-df-generate:hover {\n",
              "        background-color: #434B5C;\n",
              "        box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "        filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "        fill: #FFFFFF;\n",
              "      }\n",
              "    </style>\n",
              "    <button class=\"colab-df-generate\" onclick=\"generateWithVariable('df')\"\n",
              "            title=\"Generate code using this dataframe.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M7,19H8.4L18.45,9,17,7.55,7,17.6ZM5,21V16.75L18.45,3.32a2,2,0,0,1,2.83,0l1.4,1.43a1.91,1.91,0,0,1,.58,1.4,1.91,1.91,0,0,1-.58,1.4L9.25,21ZM18.45,9,17,7.55Zm-12,3A5.31,5.31,0,0,0,4.9,8.1,5.31,5.31,0,0,0,1,6.5,5.31,5.31,0,0,0,4.9,4.9,5.31,5.31,0,0,0,6.5,1,5.31,5.31,0,0,0,8.1,4.9,5.31,5.31,0,0,0,12,6.5,5.46,5.46,0,0,0,6.5,12Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "    <script>\n",
              "      (() => {\n",
              "      const buttonEl =\n",
              "        document.querySelector('#id_f05dd4cf-2dc8-426b-8b2b-c1d115c5d749 button.colab-df-generate');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      buttonEl.onclick = () => {\n",
              "        google.colab.notebook.generateWithVariable('df');\n",
              "      }\n",
              "      })();\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "    </div>\n",
              "  </div>\n"
            ],
            "text/plain": [
              "       p[1]    p[2]   n[1]\n",
              "0    356.12  197.67  108.0\n",
              "1    358.05  189.68   66.0\n",
              "2    340.79  260.35  130.0\n",
              "3    353.76  133.53   55.0\n",
              "4    341.37  229.80   91.0\n",
              "..      ...     ...    ...\n",
              "995  357.63  241.54   68.0\n",
              "996  352.58  212.95   87.0\n",
              "997  355.28  189.50   94.0\n",
              "998  369.75  166.33   51.0\n",
              "999  349.31  222.07  114.0\n",
              "\n",
              "[1000 rows x 3 columns]"
            ]
          },
          "execution_count": 3,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "df = pd.read_csv('https://raw.githubusercontent.com/Gurobi/modeling-examples/master/pricing_competing_products/price_value_data.csv')\n",
        "df"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b65707ac",
      "metadata": {
        "id": "b65707ac"
      },
      "source": [
        "### What's in the data?\n",
        "The data contains three columns:\n",
        "1. `p[1]` is the price (in dollars) of the first category (let's call it Category 1).\n",
        "2. `p[2]` is the price (in dollars) of the second category (Category 2).\n",
        "3. `n[1]` is the number of the items sold that are of Category 1.\n",
        "\n",
        "We don't see a column for `n[2]`, which would be the number of items sold that are Category 2. Here is where we make a **pretty big assumption** that we will sell all of the items. This makes our decision to be how to divvy up the limited space we have in order to maximize our revenue.\n",
        "The data was created to have a couple of key characteristics.\n",
        "1. As the price of Category 1 goes up, the number sold should decrease, so `p[1]` and `n[1]` have a negative correlation.\n",
        "2. As the price of Category 2 goes up, the number sold of Category 1 should increase, so `p[2]` and `n[1]` have a positive correlation.\n",
        "\n",
        "The correlation plot of the columns of the data is below."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "4cb22fe7",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 468
        },
        "id": "4cb22fe7",
        "outputId": "37183ae6-b780-4225-cc3d-bbaa2f4c7924"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAABE8AAAHDCAYAAADV6DQTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSKklEQVR4nO3de1zUVf7H8feAOoDITW7iBbyUaaYopmK1WpJ4S21t09JQ8lpZJlaGmXjZYk0z7WrmvdWfdrVNyzTTrqSmS95NTSVLQFFABQeE+f3RNruTMDLMIDK8no/H97F9z/ecM5/vPPYx5cfzOcdgNpvNAgAAAAAAQIncKjsAAAAAAACAaxnJEwAAAAAAABtIngAAAAAAANhA8gQAAAAAAMAGkicAAAAAAAA2kDwBAAAAAACwgeQJAAAAAACADSRPAAAAAAAAbCB5AgAAAAAAYAPJEwAAqpClS5fKYDDo2LFjTpvz2LFjMhgMWrp0qdPmBAAAcCUkTwAAkHTkyBGNHj1aTZo0kYeHh3x8fHTLLbdo3rx5ys/Pr+zwnGLlypWaO3duZYcBAABQ5dSo7AAAAKhs69at09/+9jcZjUbFxcWpVatWKigo0DfffKMnn3xSe/fu1YIFCyo7TIetXLlSe/bs0eOPP27VHh4ervz8fNWsWbNyAgMAALjGkTwBAFRrR48e1aBBgxQeHq4vvvhC9erVszx75JFHdPjwYa1bt86hzzCbzbp48aI8PT0ve3bx4kXVqlVLbm6VtxjUYDDIw8Oj0j4fAADgWkfZDgCgWnvhhRd0/vx5LVq0yCpx8odmzZpp3LhxkqRLly5pxowZatq0qYxGoyIiIjRp0iSZTCarMREREerTp48+++wztW/fXp6ennrzzTe1ZcsWGQwGrVq1SpMnT1b9+vXl5eWl3NxcSdLWrVvVo0cP+fr6ysvLS126dNG33357xXf46KOP1Lt3b4WFhcloNKpp06aaMWOGioqKLH26du2qdevW6fjx4zIYDDIYDIqIiJBU+p4nX3zxhW677TbVrl1bfn5+6tevn/bv32/VZ+rUqTIYDDp8+LCGDRsmPz8/+fr6Kj4+Xnl5eVZ9N27cqFtvvVV+fn7y9vZW8+bNNWnSpCu+HwAAQGVj5QkAoFr7+OOP1aRJE3Xu3PmKfUeMGKFly5bpnnvu0YQJE7R161YlJydr//79+vDDD636Hjx4UPfdd59Gjx6tkSNHqnnz5pZnM2bMUK1atfTEE0/IZDKpVq1a+uKLL9SzZ09FRUUpKSlJbm5uWrJkie644w59/fXX6tChQ6lxLV26VN7e3kpISJC3t7e++OILTZkyRbm5uZo1a5Yk6ZlnnlFOTo5OnDihl156SZLk7e1d6pyff/65evbsqSZNmmjq1KnKz8/XK6+8oltuuUU7d+60JF7+cO+996px48ZKTk7Wzp07tXDhQgUHB2vmzJmSpL1796pPnz5q3bq1pk+fLqPRqMOHD5cpOQQAAFDpzAAAVFM5OTlmSeZ+/fpdsW9qaqpZknnEiBFW7U888YRZkvmLL76wtIWHh5slmdevX2/Vd/PmzWZJ5iZNmpjz8vIs7cXFxebrrrvOHBsbay4uLra05+XlmRs3bmy+8847LW1LliwxSzIfPXrUqt+fjR492uzl5WW+ePGipa13797m8PDwy/oePXrULMm8ZMkSS1tkZKQ5ODjYnJWVZWn78ccfzW5ubua4uDhLW1JSklmS+cEHH7Sa8+677zbXrVvXcv/SSy+ZJZlPnTp12ecDAABc6yjbAQBUW3+Uy9SpU+eKfT/55BNJUkJCglX7hAkTJOmyfVEaN26s2NjYEucaOnSo1f4nqampOnTokO6//35lZWXp9OnTOn36tC5cuKBu3brpq6++UnFxcamx/e9c586d0+nTp3XbbbcpLy9PBw4cuOK7/dnJkyeVmpqqYcOGKSAgwNLeunVr3XnnnZbv4n+NGTPG6v62225TVlaW5Tv28/OT9HuJka13AQAAuBaRPAEAVFs+Pj6Sfk84XMnx48fl5uamZs2aWbWHhobKz89Px48ft2pv3LhxqXP9+dmhQ4ck/Z5UCQoKsroWLlwok8mknJycUufbu3ev7r77bvn6+srHx0dBQUEaMmSIJNkcV5o/3uV/S43+0KJFC0ti5381atTI6t7f31+SdPbsWUnSwIEDdcstt2jEiBEKCQnRoEGD9M4775BIAQAAVQJ7ngAAqi0fHx+FhYVpz549ZR5jMBjK1K+kk3VKe/ZHAmHWrFmKjIwscUxp+5NkZ2erS5cu8vHx0fTp09W0aVN5eHho586dmjhx4lVLTri7u5fYbjabJf3+zl999ZU2b96sdevWaf369Vq9erXuuOMObdiwodTxAAAA1wKSJwCAaq1Pnz5asGCBUlJSFB0dXWq/8PBwFRcX69ChQ2rRooWlPSMjQ9nZ2QoPDy93DE2bNpX0ezInJibGrrFbtmxRVlaWPvjgA/3lL3+xtB89evSyvmVN/PzxLgcPHrzs2YEDBxQYGKjatWvbFackubm5qVu3burWrZvmzJmj559/Xs8884w2b95s93sDAABcTZTtAACqtaeeekq1a9fWiBEjlJGRcdnzI0eOaN68eerVq5ckae7cuVbP58yZI0nq3bt3uWOIiopS06ZNNXv2bJ0/f/6y56dOnSp17B8rNv5Y4SFJBQUFev311y/rW7t27TKV8dSrV0+RkZFatmyZsrOzLe179uzRhg0bLN+FPc6cOXNZ2x+rbP581DMAAMC1hpUnAIBqrWnTplq5cqUGDhyoFi1aKC4uTq1atVJBQYG+++47vfvuuxo2bJjGjRunoUOHasGCBZZSmW3btmnZsmXq37+/br/99nLH4ObmpoULF6pnz5668cYbFR8fr/r16+vXX3/V5s2b5ePjo48//rjEsZ07d5a/v7+GDh2qxx57TAaDQW+//bZVMuUPUVFRWr16tRISEnTzzTfL29tbd911V4nzzpo1Sz179lR0dLSGDx9uOarY19dXU6dOtfsdp0+frq+++kq9e/dWeHi4MjMz9frrr6tBgwa69dZb7Z4PAADgaiJ5AgCo9vr27atdu3Zp1qxZ+uijj/TGG2/IaDSqdevWevHFFzVy5EhJ0sKFC9WkSRMtXbpUH374oUJDQ5WYmKikpCSHY+jatatSUlI0Y8YMvfrqqzp//rxCQ0PVsWNHjR49utRxdevW1dq1azVhwgRNnjxZ/v7+GjJkiLp163bZaT8PP/ywUlNTtWTJEr300ksKDw8vNXkSExOj9evXKykpSVOmTFHNmjXVpUsXzZw50+ZmuKXp27evjh07psWLF+v06dMKDAxUly5dNG3aNPn6+to9HwAAwNVkMJf0V1MAAAAAAACQxJ4nAAAAAAAANpE8AQAAAAAAsIHkCQAAAAAAgA0kTwAAAAAAgNN99dVXuuuuuxQWFiaDwaA1a9ZcccyWLVvUrl07GY1GNWvWTEuXLr2sz2uvvaaIiAh5eHioY8eO2rZtm/OD/xOSJwAAAAAAwOkuXLigNm3a6LXXXitT/6NHj6p37966/fbblZqaqscff1wjRozQZ599ZumzevVqJSQkKCkpSTt37lSbNm0UGxurzMzMinoNSZy2AwAAAAAAKpjBYNCHH36o/v37l9pn4sSJWrdunfbs2WNpGzRokLKzs7V+/XpJUseOHXXzzTfr1VdflSQVFxerYcOGevTRR/X0009XWPysPAEAAAAAAGViMpmUm5trdZlMJqfMnZKSopiYGKu22NhYpaSkSJIKCgq0Y8cOqz5ubm6KiYmx9KkoNSp0dnvsfb+yIwAAp4ro9URlhwAATvXYqjWVHQIAOFVCdJvKDuHqcOKft5Pf3a1p06ZZtSUlJWnq1KkOz52enq6QkBCrtpCQEOXm5io/P19nz55VUVFRiX0OHDjg8Ofbcu0kTwAAAAAAwDUtMTFRCQkJVm1Go7GSorl6SJ4AAAAAAIAyMRqNFZYsCQ0NVUZGhlVbRkaGfHx85OnpKXd3d7m7u5fYJzQ0tEJi+gN7ngAAAAAA4MLMRUVOuypSdHS0Nm3aZNW2ceNGRUdHS5Jq1aqlqKgoqz7FxcXatGmTpU9FIXkCAAAAAACc7vz580pNTVVqaqqk348iTk1NVVpamqTfS4Di4uIs/ceMGaOff/5ZTz31lA4cOKDXX39d77zzjsaPH2/pk5CQoLfeekvLli3T/v379dBDD+nChQuKj4+v0HehbAcAAAAAAFdWdKlSPvaHH37Q7bffbrn/Y6+UoUOHaunSpTp58qQlkSJJjRs31rp16zR+/HjNmzdPDRo00MKFCxUbG2vpM3DgQJ06dUpTpkxRenq6IiMjtX79+ss2kXU2g9lsNlfoJ5QVp+0AcDGctgPA1XDaDgBXU11O2ynesdxpc7lFxV25kwuibAcAAAAAAMAGynYAAAAAAHBlFbzRa3VA8gQAAAAAABdmrqQ9T1wJZTsAAAAAAAA2sPIEAAAAAABXxsoTh5E8AQAAAADAhZmLSZ44irIdAAAAAAAAG1h5AgAAAACAK+O0HYeRPAEAAAAAwIVx2o7jKNsBAAAAAACwgZUnAAAAAAC4MlaeOIzkCQAAAAAALsxczJ4njqJsBwAAAAAAwAZWngAAAAAA4MLYMNZxJE8AAAAAAHBlJE8cRtkOAAAAAACADaw8AQAAAADAhbFhrONYeQIAAAAAAGADyRMAAAAAAAAbKNsBAAAAAMCVsWGsw0ieAAAAAADgwjiq2HGU7QAAAAAAANjAyhMAAAAAAFwZK08cRvIEAAAAAAAXxlHFjqNsBwAAAAAAwAZWngAAAAAA4Moo23EYyRMAAAAAAFyYuYiyHUdRtgMAAAAAAGADK08AAAAAAHBhZsp2HEbyBAAAAAAAV1ZM8sRRlO0AAAAAAADYwMoTAAAAAABcGBvGOo7kCQAAAAAArozkicMo2wEAAAAAALCBlScAAAAAALgwTttxHMkTAAAAAABcGWU7DqNsBwAAAAAAwAZWngAAAAAA4MI4bcdxJE8AAAAAAHBh5mKSJ46ibAcAAAAAAMAGkicAAAAAALiyoiLnXeXw2muvKSIiQh4eHurYsaO2bdtWat+uXbvKYDBcdvXu3dvSZ9iwYZc979GjR7liKyvKdgAAAAAAQIVYvXq1EhISNH/+fHXs2FFz585VbGysDh48qODg4Mv6f/DBByooKLDcZ2VlqU2bNvrb3/5m1a9Hjx5asmSJ5d5oNFbcS4iVJwAAAAAAoILMmTNHI0eOVHx8vFq2bKn58+fLy8tLixcvLrF/QECAQkNDLdfGjRvl5eV1WfLEaDRa9fP396/Q9yB5AgAAAACACzMXFTntskdBQYF27NihmJgYS5ubm5tiYmKUkpJSpjkWLVqkQYMGqXbt2lbtW7ZsUXBwsJo3b66HHnpIWVlZdsVmL8p2AAAAAABwYeaiYqfNZTKZZDKZrNqMRmOJZTOnT59WUVGRQkJCrNpDQkJ04MCBK37Wtm3btGfPHi1atMiqvUePHvrrX/+qxo0b68iRI5o0aZJ69uyplJQUubu7l+OtroyVJwAAAAAAoEySk5Pl6+trdSUnJ1fIZy1atEg33XSTOnToYNU+aNAg9e3bVzfddJP69++vtWvXavv27dqyZUuFxCGRPAEAAAAAwLUVFTvtSkxMVE5OjtWVmJhY4scGBgbK3d1dGRkZVu0ZGRkKDQ21GfKFCxe0atUqDR8+/Iqv16RJEwUGBurw4cNl/07sVOaynYSEBLsnnzx5sgICAuweBwAAAAAAnMPevUpsKa1EpyS1atVSVFSUNm3apP79+0uSiouLtWnTJo0dO9bm2HfffVcmk0lDhgy54uecOHFCWVlZqlevXpniKo8yJ0/mzp2r6Oho1apVq0z9v/nmG40dO5bkCQAAAAAA1VRCQoKGDh2q9u3bq0OHDpo7d64uXLig+Ph4SVJcXJzq169/WenPokWL1L9/f9WtW9eq/fz585o2bZoGDBig0NBQHTlyRE899ZSaNWum2NjYCnsPuzaM/fDDD0s8h7kkderUKVdAAAAAAADAecxF5kr77IEDB+rUqVOaMmWK0tPTFRkZqfXr11s2kU1LS5Obm/WOIgcPHtQ333yjDRs2XDafu7u7du3apWXLlik7O1thYWHq3r27ZsyYUeYVMeVR5uTJkiVL5OvrW+aJ33zzzct21AUAAAAAAFeXM0/bKY+xY8eWWqZT0iavzZs3l9lccsLH09NTn332mTPDK5MyJ0+GDh1q18T333+/3cEAFWX73qNa9NHX2nPkV506e06vTRyimI4tKzssACjV+ITxuu++QfLx8dEPP/ygyc88q2PHjtkc80DcAxo9apSCgoK0f/9+JSVN1Y8//ihJ8vX11fiE8brttttUv36YsrKytGHDRs15cY7OnTt3Fd4IQHVlNpv1w4fv6MCXm2TKu6DQ627QbXEj5Bta+t4ExcXF2vHhOzqU8rXycrJV2y9A19/aRe36DpDBYLD0O/vbCW19Z4VOHtyn4qJi+ddvoDvHTlCduoFX49UAVCOctoNqIc9UoOYRoUoa2beyQwGAKxozZrTihw3TM5Mmq3+/u5Wfl6/lby+T0Vj6vmN9+vTW5MnPaN68eerdp4/27d+v5W8vs9QJh4SEKCQkWM8/97y63xmrJ554Ul26dNHMF2ZerdcCUE39+MlH2rPxU902dKTunvK8ahiNWvfic7pUUFDqmNR1a7Rv80bdMmS4Bj7/kjreO1g/fvov7fn8U0ufnMx0ffTcFPnVq6+7np6qe/4+S+36DlCNmjWvxmsBVYq5qNhpV3Xl1OTJjz/+KHd3d2dOCThFl3bNNf7+7rqz042VHQoAXNGDwx/UK6++qo0bN+rAgQNKSJigkOAQde/evdQxI0aM0KpVq/Xuu+/p8KHDembSM8rPz9e99/5NkvTTTz/poTEPa9OmTUpLS1PKdymaPWu2unW7g393A6gwZrNZuzd8onZ9/6qIdjerbsNw3T5yrPLOntWxndtLHZdx+CeFt22v8Mh2qhMUrCY3d1KDG1sr8+f/HkO6/b1VatS6rToNHKLA8MbyDQ5VRNv28vQp+1YDQHVhLjY77aqunL7ypLS6JAAAcGUNGzZUcHCwvv3mG0vbuXPnlJqaqnbt2pU4pmbNmmp1UyurMWazWd9+822pYySpjk8dnT9/XkVOPL4QAP7XuVOZysvJVv2WrS1tRi8vBTdtpowjP5U6LqTZ9fp13x5lp/8mScpKO6b0QwfV6Ka2kiRzcbHSdu2Ub2g9rZv9nJY9OkIfTp+kozu2VewLAai27Dpt569//avN5zk5OVY1iAAAwD5BwUGSpFOnT1u1nzp9WkFBQSWO8ff3V40aNXS6hDFNmzYtdcyjjz6q//u/VU6IGgBKlpeTLUny/NPBE54+vpZnJWnbu78K8/O1OnG83NzcVFxcrA4DBum6zrdJkvJzc1V48aJS132kmwcMVMe/DdYvu1O14dUXddfEJIXdwN52wP+qzNN2XIVdyZOPP/5Yd955Z6mn6JT1b65MJpNMJpNVm7GgUMZa1CcCAKqXfv376fnnn7PcPxg/vMI/09vbW0uWLNbhw4c096W5Ff55AKqPQ999ra+WLbDc9xyfWK55jmxL0aHvv1G30Y/Jv35DZaUd03crl8rLz1/Nb+0qs/n3fRci2rVX69g+kqTA8AhlHD6ofZs3kDwB/sTMIlOH2ZU8adGihQYMGKDhw0v+D7vU1FStXbv2ivMkJydr2rRpVm1JD/1NUx8ZaE84AABUeZ9v/Fyp/0613Neq9fumsEGBgTqVecrSHhQYqH379pU4x9mzZ3Xp0iUFBlqfLhEUGKhTp05ZtdWuXVvLli/V+QvnNXrUaF26dMlJbwIAUnjb9rqn6XWW+6JLhZKk/Jwc1fbzt7Tn5+aobqOIUuf5/p1/KrJXPzXrdIskqW7DRjqfdUqpa9eo+a1d5VHHR27u7vIPa2A1zi+svtJ/OujENwKA39m150lUVJR27txZ6nOj0ahGjRpdcZ7ExETl5ORYXYkjbZcEAQDgii5cuKDjx49brkOHDikzM1Odb7nF0sfb21uRkZGl/ju4sLBQe3bvsRpjMBjU+ZbOVmO8vb319j+Xq7CgUCOGj5TJVPpJFwBQHrU8PeUbEmq5/MMayMvXT7/u223pU5Cfp8wjhxXS9PpS57lkMsngZv1HFYObm2V/RfcaNRTUuKmyT/5m1Scn/aTqBHJMMfBn5iKz067qyq6VJ/Pnz7dZmtOiRQsdPXr0ivMYjUYZjUbrRkp2UIEu5JuUlp5luT+ReUb7j/4mX28vhQX5VV5gAFCCxYsW69FHx+rY0WP65ZdfNGFCgjIyM7RhwwZLnxUr/6nPPtug5cuWS5IWLlyoF198Ubt37VLqjz9q+IMPysvLS++++56k/yRO3l4uD09PPT5uvOrU8VadOt6SpKysMyourr5HDwKoOAaDQTd176WdH38g39B6qhMYrB8+WCUvf39FtLvZ0u/jmdPVOKqDWsX0kCSFR0bp3x9/IO+AQAXUb6DTace067O1an7b7ZYxbXr21eevv6R6zVsorEUr/bI7VcdTd+iup6de7dcErnn8a95xdiVPLkt4AFXEniO/Km7KQst98pJPJEl3395O/3j0nsoKCwBKNH/+m/L08lJy8vPy8fHR9h+2a2jcMKuVIuGNwhXg/98l8GvXrlNA3boan5CgoKBA7d+3X0Pjhlk2kW3V6ka1bff7KRVfff2l1efdesutOnHi16vwZgCqoza9+qnQZNJXS95UQV6eQq+/Qb0mTFKN/5QpSlJuZoYunsu13N8y5EFt/2C1vnl7ofJzc1TbL0Atut6pqH7//e+2xlEddNvQkfr3ujX6dsUS+YWGqfvYCap3/Q1X9f0AVA8GcxnPFs7NzZWPj0+ZJz537pzq1KlT9kj2vl/2vgBQBUT0eqKyQwAAp3ps1ZrKDgEAnCohuk1lh3BV/DKsvdPmarj0B6fNVZWUec8Tf39/ZWZmlnni+vXr6+effy5XUAAAAAAAwDnMRc67qqsyl+2YzWYtXLhQ3t7eZepfWFhY7qAAAAAAAACuFWVOnjRq1EhvvfVWmScODQ1VzZpsAgsAAAAAAKq2MidPjh07VmL7H1umGAwGpwQEAAAAAACch9N2HFfmPU/+bNGiRWrVqpU8PDzk4eGhVq1aaeHChVceCAAAAAAAUIXYdVTxH6ZMmaI5c+bo0UcfVXR0tCQpJSVF48ePV1pamqZPn+7UIAEAAAAAQPlU541enaVcyZM33nhDb731lu677z5LW9++fdW6dWs9+uijJE8AAAAAALhGFBezzYajylW2U1hYqPbtLz8nOioqSpcuXXI4KAAAAAAAgGtFuZInDzzwgN54443L2hcsWKDBgwc7HBQAAAAAAHCO4mLnXdVVucp2pN83jN2wYYM6deokSdq6davS0tIUFxenhIQES785c+Y4HiUAAAAAACgX9jxxXLmSJ3v27FG7du0kSUeOHJEkBQYGKjAwUHv27LH04/hiAAAAAABQ1ZUrebJ582ZnxwEAAAAAACoAG8Y6rtxlOwAAAAAA4NpXTNmOw8q1YSwAAAAAAEB1wcoTAAAAAABcGGU7jiN5AgAAAACACzOTPHEYZTsAAAAAAAA2sPIEAAAAAAAXVlxc2RFUfaw8AQAAAAAAsIGVJwAAAAAAuDA2jHUcyRMAAAAAAFwYyRPHUbYDAAAAAABgAytPAAAAAABwYUWsPHEYyRMAAAAAAFwYZTuOo2wHAAAAAADABlaeAAAAAADgworNrDxxFCtPAAAAAAAAbCB5AgAAAAAAYANlOwAAAAAAuLDi4sqOoOojeQIAAAAAgAsrYs8Th1G2AwAAAAAAKsxrr72miIgIeXh4qGPHjtq2bVupfZcuXSqDwWB1eXh4WPUxm82aMmWK6tWrJ09PT8XExOjQoUMV+g4kTwAAAAAAcGHFxQanXfZavXq1EhISlJSUpJ07d6pNmzaKjY1VZmZmqWN8fHx08uRJy3X8+HGr5y+88IJefvllzZ8/X1u3blXt2rUVGxurixcv2h1fWZE8AQAAAADAhRWZDU677DVnzhyNHDlS8fHxatmypebPny8vLy8tXry41DEGg0GhoaGWKyQkxPLMbDZr7ty5mjx5svr166fWrVtr+fLl+u2337RmzZryfD1lQvIEAAAAAAA4XUFBgXbs2KGYmBhLm5ubm2JiYpSSklLquPPnzys8PFwNGzZUv379tHfvXsuzo0ePKj093WpOX19fdezY0eacjmLDWAAAAAAAXFixEzeMNZlMMplMVm1Go1FGo/GyvqdPn1ZRUZHVyhFJCgkJ0YEDB0qcv3nz5lq8eLFat26tnJwczZ49W507d9bevXvVoEEDpaenW+b485x/PKsIrDwBAAAAAMCFObNsJzk5Wb6+vlZXcnKy02KNjo5WXFycIiMj1aVLF33wwQcKCgrSm2++6bTPKA9WngAAAAAAgDJJTExUQkKCVVtJq04kKTAwUO7u7srIyLBqz8jIUGhoaJk+r2bNmmrbtq0OHz4sSZZxGRkZqlevntWckZGRZX0Nu7HyBAAAAAAAF1Zkdt5lNBrl4+NjdZWWPKlVq5aioqK0adMmS1txcbE2bdqk6OjossVeVKTdu3dbEiWNGzdWaGio1Zy5ubnaunVrmecsD1aeAAAAAADgwpy554m9EhISNHToULVv314dOnTQ3LlzdeHCBcXHx0uS4uLiVL9+fUvpz/Tp09WpUyc1a9ZM2dnZmjVrlo4fP64RI0ZI+v0knscff1x///vfdd1116lx48Z69tlnFRYWpv79+1fYe5A8AQAAAAAAFWLgwIE6deqUpkyZovT0dEVGRmr9+vWWDV/T0tLk5vbfopizZ89q5MiRSk9Pl7+/v6KiovTdd9+pZcuWlj5PPfWULly4oFGjRik7O1u33nqr1q9fLw8Pjwp7D4PZbDZX2Oz22Pt+ZUcAAE4V0euJyg4BAJzqsVVrKjsEAHCqhOg2lR3CVfFhu5grdyqju3d+7rS5qhJWngAAAAAA4MKKro0lE1UaG8YCAAAAAADYwMoTAAAAAABcWJEqb8NYV0HyBAAAAAAAF0bZjuMo2wEAAAAAALCBlScAAAAAALiwosoOwAWQPAEAAAAAwIWRPHEcZTsAAAAAAAA2kDwBAAAAAACwgbIdAAAAAABcGEcVO46VJwAAAAAAADaw8gQAAAAAABdWZDZXdghVHskTAAAAAABcGKftOI6yHQAAAAAAABtYeQIAAAAAgAtj5YnjSJ4AAAAAAODCSJ44jrIdAAAAAAAAG1h5AgAAAACACysSp+04iuQJAAAAAAAujLIdx10zyZOIXk9UdggA4FTHPpld2SEAgFP17D+kskMAAKdKOLS7skNAFXHNJE8AAAAAAIDzFZkp23EUyRMAAAAAAFwYZTuO47QdAAAAAAAAG1h5AgAAAACAC+O0HceRPAEAAAAAwIWRPHEcZTsAAAAAAAA2sPIEAAAAAAAXxoaxjiN5AgAAAACAC+OoYsdRtgMAAAAAAGADK08AAAAAAHBhbBjrOJInAAAAAAC4MJInjqNsBwAAAAAAwAaSJwAAAAAAADZQtgMAAAAAgAsr5rQdh7HyBAAAAAAAwAZWngAAAAAA4MLYMNZxJE8AAAAAAHBhJE8cR9kOAAAAAACADaw8AQAAAADAhRWxYazDSJ4AAAAAAODCKNtxHGU7AAAAAAAANpA8AQAAAADAhRWbzU67yuO1115TRESEPDw81LFjR23btq3Uvm+99ZZuu+02+fv7y9/fXzExMZf1HzZsmAwGg9XVo0ePcsVWViRPAAAAAABwYUUyO+2y1+rVq5WQkKCkpCTt3LlTbdq0UWxsrDIzM0vsv2XLFt13333avHmzUlJS1LBhQ3Xv3l2//vqrVb8ePXro5MmTluv//u//yvXdlBXJEwAAAAAAUCHmzJmjkSNHKj4+Xi1bttT8+fPl5eWlxYsXl9h/xYoVevjhhxUZGakbbrhBCxcuVHFxsTZt2mTVz2g0KjQ01HL5+/tX6HuQPAEAAAAAwIU5c+WJyWRSbm6u1WUymUr83IKCAu3YsUMxMTGWNjc3N8XExCglJaVMsefl5amwsFABAQFW7Vu2bFFwcLCaN2+uhx56SFlZWeX/gsqA5AkAAAAAAC7MmXueJCcny9fX1+pKTk4u8XNPnz6toqIihYSEWLWHhIQoPT29TLFPnDhRYWFhVgmYHj16aPny5dq0aZNmzpypL7/8Uj179lRRUVH5v6Qr4KhiAAAAAABQJomJiUpISLBqMxqNFfJZ//jHP7Rq1Spt2bJFHh4elvZBgwZZ/vmmm25S69at1bRpU23ZskXdunWrkFhIngAAAAAA4MLKs9FraYxGY5mTJYGBgXJ3d1dGRoZVe0ZGhkJDQ22OnT17tv7xj3/o888/V+vWrW32bdKkiQIDA3X48OEKS55QtgMAAAAAgAsrMpuddtmjVq1aioqKstrs9Y/NX6Ojo0sd98ILL2jGjBlav3692rdvf8XPOXHihLKyslSvXj274rMHyRMAAAAAAFAhEhIS9NZbb2nZsmXav3+/HnroIV24cEHx8fGSpLi4OCUmJlr6z5w5U88++6wWL16siIgIpaenKz09XefPn5cknT9/Xk8++aS+//57HTt2TJs2bVK/fv3UrFkzxcbGVth7ULYDAAAAAIALK3Zi2Y69Bg4cqFOnTmnKlClKT09XZGSk1q9fb9lENi0tTW5u/13X8cYbb6igoED33HOP1TxJSUmaOnWq3N3dtWvXLi1btkzZ2dkKCwtT9+7dNWPGjArbe0WSDGaznetuKkhEeOPKDgEAnOrYJ7MrOwQAcKqe/adWdggA4FSfHtpd2SFcFf2bt3XaXGsO/ttpc1UllO0AAAAAAADYQNkOAAAAAAAurPjaKDip0lh5AgAAAAAAYAPJEwAAAAAAABso2wEAAAAAwIUVVeJpO66C5AkAAAAAAC6s2Fxc2SFUeZTtAAAAAAAA2MDKEwAAAAAAXFgxZTsOI3kCAAAAAIALK+KoYodRtgMAAAAAAGADK08AAAAAAHBhlO04juQJAAAAAAAurJiyHYdRtgMAAAAAAGADK08AAAAAAHBhxZUdgAsgeQIAAAAAgAujbMdxlO0AAAAAAADYwMoTAAAAAABcGKftOI7kCQAAAAAALoyyHcdRtgMAAAAAAGADK08AAAAAAHBhlO04rszJk5dfftnuyePj41WnTh27xwEAAAAAAOcgeeK4MidPHn/8cTVo0EDu7u5l6v/LL7+oT58+JE8AAAAAAECVZlfZzg8//KDg4OAy9SVpAgAAAABA5Stm4YnDypw8SUpKkre3d5knnjRpkgICAsoVFAAAAAAAcA7KdhxnV/LEHomJiXYHAwAAAAAAcK3hqGK4hPEJ47Vt+1YdOLhf/1zxtiIiIq445oG4B/TNN1/r4MEDWrPmQ7Vp08byzNfXV1OnTdWmLzbpwMH9+va7b5Q0NYlyNADXjO17j2rM88t16/BkNf/rJH2+dV9lhwQAJeozeJCWbl6vj/b8oJfeW6HrW7cqtW+jZk31zKtztHTzen16aLf6DxtyWR/P2l4a/cxTWrrlM63ZvV0vrn5b1990Y0W+AlDlFcvstKu6sjt58sknn2jEiBF66qmndODAAatnZ8+e1R133OG04ICyGDNmtOKHDdMzkyarf7+7lZ+Xr+VvL5PRWKvUMX369Nbkyc9o3rx56t2nj/bt36/lby9T3bp1JUkhISEKCQnW8889r+53xuqJJ55Uly5dNPOFmVfrtQDApjxTgZpHhCppZN/KDgUASvWXXrEaNelJrXh1vh7tf6+O7v9Jf1/8pnxLKe/38PRQ+i8ntGT2XJ3JPFVin3HPTVPbW6I1+8lJeqj3X7Xzm+/0/LK3VDekbHszAkB52JU8Wblypfr27av09HSlpKSobdu2WrFiheV5QUGBvvzyS6cHCdjy4PAH9cqrr2rjxo06cOCAEhImKCQ4RN27dy91zIgRI7Rq1Wq9++57OnzosJ6Z9Izy8/N1771/kyT99NNPemjMw9q0aZPS0tKU8l2KZs+arW7d7ijziVMAUJG6tGuu8fd3152d+NtWANeuux+M06er39fG99co7fDPemXKdJny89X9nrtL7P/T7r1aNHOOvly3XoUFBZc9r2U06tbYGC16YY72bN+hk2m/aMUrb+i347+o9/0DK/p1AFRjdiVPZs2apTlz5mjt2rX6+uuvtWzZMo0ePVqLFi2qqPgAmxo2bKjg4GB9+803lrZz584pNTVV7dq1K3FMzZo11eqmVlZjzGazvv3m21LHSFIdnzo6f/68ioqKnPcCAAAALqpGzRq67saWSv3ue0ub2WxW6nffq0XbNjZGls69hrvca9RQock6sVJw8aJujGrrULyAKzObnXdVV3YdVXzo0CHdddddlvt7771XQUFB6tu3rwoLC3X33SVnkIGKEhQcJEk6dfq0Vfup06cVFBRU4hh/f3/VqFFDp0sY07Rp01LHPProo/q//1vlhKgBAABcn4+/v9xr1NDZ01lW7WezstSgaeNyzZl/IU/7dqbqvkdGK+3Iz8o+naUufXrphrZtdPJ4mjPCBlxSdd6rxFnsWnni4+OjjIwMq7bbb79da9eu1ZNPPqlXXnmlTPOYTCbl5uZaXebqnMJCmfXr30979+2xXDVr1Kzwz/T29taSJYt1+PAhzX1pboV/HgAAAEo3+8lEGQwGrfj2C/1r7w71i7tfX679VMX8eQJABbJr5UmHDh306aefqlOnTlbtXbp00ccff6w+ffqUaZ7k5GRNmzbNqs3Xx1d+fv72hINq6PONnyv136mW+1q1ft8UNigwUKf+Z1OxoMBA7dtX8skTZ8+e1aVLlxQYGGjVHhQYqFOnrDcmq127tpYtX6rzF85r9KjRunTpkpPeBAAAwLXlnj2rokuX5B9Y16rdv25dnT2VVcqoKzuZdkJPDY6X0dNTXt61dfbUaT09d5bSfznhaMiAyyK16Di7Vp6MHz9eHh4eJT7r2rWrPv74Y8XFxV1xnsTEROXk5Fhdvr5+9oSCaurChQs6fvy45Tp06JAyMzPV+ZZbLH28vb0VGRmpnTt3ljhHYWGh9uzeYzXGYDCo8y2drcZ4e3vr7X8uV2FBoUYMHymT6fJNywAAAFCyS4WXdGjvPkVGd7S0GQwGRXbupP3//tHh+U35+Tp76rS8fXwUdVtnff/5ZofnBFwVRxU7zq6VJ126dFGXLl1KfX777bfr9ttvv+I8RqNRRqPRqs1gMNgTCmCxeNFiPfroWB07eky//PKLJkxIUEZmhjZs2GDps2LlP/XZZxu0fNlySdLChQv14osvaveuXUr98UcNf/BBeXl56d1335P0n8TJ28vl4empx8eNV5063qpTx1uSlJV1RsXFxVf/RQHgf1zINykt/b9/c3si84z2H/1Nvt5eCgvyq7zAAOB/fLh4uSa88JwO7dmrg7t2q/+wB2T09NTG99dIkia88JyyMjK19MV5kn7fZLZRs6b/+eeaqhsSrCYtmiv/Qp5Opv0iSWp3a2cZDAadOHpMYeGNNHxigk78fFQb/jMnAFSEMidPcnNz5ePjU+aJz507pzp16pQrKMAe8+e/KU8vLyUnPy8fHx9t/2G7hsYNs1opEt4oXAH+/y0LW7t2nQLq1tX4hAQFBQVq/779Gho3zLKJbKtWN6ptu993bP/qa+vjt2+95VadOPHrVXgzACjdniO/Km7KQst98pJPJEl3395O/3j0nsoKCwCsfPXJZ/INCNCQcY8oIChQR/Yf0LPDxyg76/fkb3BYPau9DwOCg/Xav96z3N8zIl73jIjXrq3bNXHIg5Kk2nXqKP6JcQoMDdG57Bx989nnWjbnZRVRXg2UqvquF3Eeg7mMO7W6u7vr5MmTCg4OLtPEPj4+Sk1NVZMmTcrUPyK8fDtuA8C16tgnsys7BABwqp79p1Z2CADgVJ8e2l3ZIVwV14dHOG2un44fc9pcVUmZV56YzWYtXLhQ3t7eZepfWFhY7qAAAAAAAACuFWVOnjRq1EhvvfVWmScODQ1VzZoVf4wsAAAAAAAoXXXe6NVZypw8OXbsWIntf1T9sOErAAAAAADXHlInjrPrqOL/tWjRIrVq1UoeHh7y8PBQq1attHDhwisPBAAAAAAAqELKlTyZMmWKxo0bp7vuukvvvvuu3n33Xd11110aP368pkyZ4uwYAQAAAABAOZmdeJXHa6+9poiICHl4eKhjx47atm2bzf7vvvuubrjhBnl4eOimm27SJ598Yv0+ZrOmTJmievXqydPTUzExMTp06FA5oyubciVP3njjDb311ltKTk5W37591bdvXyUnJ2vBggV6/fXXnR0jAAAAAAAop8pMnqxevVoJCQlKSkrSzp071aZNG8XGxiozM7PE/t99953uu+8+DR8+XP/+97/Vv39/9e/fX3v27LH0eeGFF/Tyyy9r/vz52rp1q2rXrq3Y2FhdvHixHBGWTbmSJ4WFhWrfvv1l7VFRUbrE+eoAAAAAAEDSnDlzNHLkSMXHx6tly5aaP3++vLy8tHjx4hL7z5s3Tz169NCTTz6pFi1aaMaMGWrXrp1effVVSb+vOpk7d64mT56sfv36qXXr1lq+fLl+++03rVmzpsLeo1zJkwceeEBvvPHGZe0LFizQ4MGDHQ4KAAAAAAA4R2WtPCkoKNCOHTsUExNjaXNzc1NMTIxSUlJKHJOSkmLVX5JiY2Mt/Y8ePar09HSrPr6+vurYsWOpczpDmU/b+bNFixZpw4YN6tSpkyRp69atSktLU1xcnBISEiz95syZ43iUAAAAAACg0plMJplMJqs2o9Eoo9F4Wd/Tp0+rqKhIISEhVu0hISE6cOBAifOnp6eX2D89Pd3y/I+20vpUhHIlT/bs2aN27dpJko4cOSJJCgwMVGBgoFUdEscXAwAAAADgOpKTkzVt2jSrtqSkJE2dOrVyArpKypU82bx5s7PjAAAAAAAAFcJ5CxsSExOtqk0klbjqRPp9kYW7u7syMjKs2jMyMhQaGlrimNDQUJv9//jfjIwM1atXz6pPZGSkXe9ij3LteQIAAAAAAKoKg9Muo9EoHx8fq6u05EmtWrUUFRWlTZs2WdqKi4u1adMmRUdHlzgmOjraqr8kbdy40dK/cePGCg0NteqTm5urrVu3ljqnM5R7zxMAAAAAAABbEhISNHToULVv314dOnTQ3LlzdeHCBcXHx0uS4uLiVL9+fSUnJ0uSxo0bpy5duujFF19U7969tWrVKv3www9asGCBpN+3B3n88cf197//Xdddd50aN26sZ599VmFhYerfv3+FvQfJEwAAAAAAUCEGDhyoU6dOacqUKUpPT1dkZKTWr19v2fA1LS1Nbm7/LYrp3LmzVq5cqcmTJ2vSpEm67rrrtGbNGrVq1crS56mnntKFCxc0atQoZWdn69Zbb9X69evl4eFRYe9hMJvN9p42VCEiwhtXdggA4FTHPpld2SEAgFP17D+1skMAAKf69NDuyg7hqogIb+K0uY4d/9lpc1Ul7HkCAAAAAABgA2U7AAAAAAC4MucdtlNtkTwBAAAAAMClUXTiKL5BAAAAAAAAG1h5AgAAAACACzNQt+MwkicAAAAAALgyA8kTR1G2AwAAAAAAYAMrTwAAAAAAcGGU7TiO5AkAAAAAAC6NohNH8Q0CAAAAAADYwMoTAAAAAABcmIENYx1G8gQAAAAAAFdmoOjEUXyDAAAAAAAANrDyBAAAAAAAF2Zg3YTDSJ4AAAAAAODC2PPEcaSfAAAAAAAAbGDlCQAAAAAArowNYx1G8gQAAAAAABdmIHniML5BAAAAAAAAG1h5AgAAAACAC+O0HceRPAEAAAAAwIVRtuM4vkEAAAAAAAAbSJ4AAAAAAADYQNkOAAAAAAAuzGBwr+wQqjxWngAAAAAAANjAyhMAAAAAAFwYG8Y6juQJAAAAAAAujOSJ4/gGAQAAAAAAbGDlCQAAAAAALowNYx1H8gQAAAAAABdG2Y7j+AYBAAAAAABsYOUJAAAAAAAujLIdx5E8AQAAAADAhZE8cRxlOwAAAAAAADaw8gQAAAAAABfmxoaxDiN5AgAAAACAC6Nsx3GknwAAAAAAAGxg5QkAAAAAAC6MlSeOI3kCAAAAAIALI3niOMp2AAAAAAAAbGDlCQAAAAAALszgxsoTR7HyBAAAAAAAF+ZmcHfaVVHOnDmjwYMHy8fHR35+fho+fLjOnz9vs/+jjz6q5s2by9PTU40aNdJjjz2mnJwcq34Gg+Gya9WqVXbHx8oTAAAAAABQqQYPHqyTJ09q48aNKiwsVHx8vEaNGqWVK1eW2P+3337Tb7/9ptmzZ6tly5Y6fvy4xowZo99++03vvfeeVd8lS5aoR48elns/Pz+74yN5AgAAAACAC7vWN4zdv3+/1q9fr+3bt6t9+/aSpFdeeUW9evXS7NmzFRYWdtmYVq1a6f3337fcN23aVM8995yGDBmiS5cuqUaN/6Y7/Pz8FBoa6lCM10zy5LFVayo7BABwqp79h1R2CADgVJ+umVrZIQAAKpnJZJLJZLJqMxqNMhqN5Z4zJSVFfn5+lsSJJMXExMjNzU1bt27V3XffXaZ5cnJy5OPjY5U4kaRHHnlEI0aMUJMmTTRmzBjFx8fLYDDYFSN7ngAAAAAAgDJJTk6Wr6+v1ZWcnOzQnOnp6QoODrZqq1GjhgICApSenl6mOU6fPq0ZM2Zo1KhRVu3Tp0/XO++8o40bN2rAgAF6+OGH9corr9gd4zWz8gQAAAAAADifM8t2EhMTlZCQYNVW2qqTp59+WjNnzrQ53/79+x2OKTc3V71791bLli01depUq2fPPvus5Z/btm2rCxcuaNasWXrsscfs+gySJwAAAAAAuDCDwXl/9LenRGfChAkaNmyYzT5NmjRRaGioMjMzrdovXbqkM2fOXHGvknPnzqlHjx6qU6eOPvzwQ9WsWdNm/44dO2rGjBkymUx2lRqRPAEAAAAAAE4XFBSkoKCgK/aLjo5Wdna2duzYoaioKEnSF198oeLiYnXs2LHUcbm5uYqNjZXRaNS//vUveXh4XPGzUlNT5e/vb/ceLSRPAAAAAABwYW7X+Gk7LVq0UI8ePTRy5EjNnz9fhYWFGjt2rAYNGmQ5aefXX39Vt27dtHz5cnXo0EG5ubnq3r278vLy9M9//lO5ubnKzc2V9HvSxt3dXR9//LEyMjLUqVMneXh4aOPGjXr++ef1xBNP2B0jyRMAAAAAAFyYwe3aTp5I0ooVKzR27Fh169ZNbm5uGjBggF5++WXL88LCQh08eFB5eXmSpJ07d2rr1q2SpGbNmlnNdfToUUVERKhmzZp67bXXNH78eJnNZjVr1kxz5szRyJEj7Y6P5AkAAAAAAKhUAQEBWrlyZanPIyIiZDabLfddu3a1ui9Jjx491KNHD6fER/IEAAAAAAAX5swNY6srvkEAAAAAAFyYM48qrq7cKjsAAAAAAACAaxkrTwAAAAAAcGGU7TiObxAAAAAAABd2rR9VXBVQtgMAAAAAAGADK08AAAAAAHBhBjf+6O8ovkEAAAAAAFwYe544jrIdAAAAAAAAG0g/AQAAAADgwgxsGOswkicAAAAAALgwynYcR9kOAAAAAACADaSfAAAAAABwYZy24zhWngAAAAAAANhA+gkAAAAAABfGnieO4xsEAAAAAMCVkTxxGGU7AAAAAAAANpA8AQAAAAAAsIG1OwAAAAAAuDBO23EcK08AAAAAAABsIP0EAAAAAIAL47Qdx/ENAgAAAADgyijbcRhlOwAAAAAAADaQfgIAAAAAwJUZ3Cs7giqP5AkAAAAAAC6M03YcR9kOAAAAAACADaSfAAAAAABwZZy24zC+QQAAAAAAXJiZsh2HUbYDAAAAAABgA+knAAAAAABcmRun7TiK5AkAAAAAAK6M5InDKNsBAAAAAACwgZUnAAAAAAC4MDMrTxxG8gQAAAAAABdG8sRxlO0AAAAAAADYwMoTAAAAAABcGStPHEbyBAAAAAAAF2Z2o+jEUXyDAAAAAAAANrDyBAAAAAAAF8aGsY5j5QkAAAAAAIANJE8AAAAAAABsIHkCAAAAAIALK3Z3c9pVUc6cOaPBgwfLx8dHfn5+Gj58uM6fP29zTNeuXWUwGKyuMWPGWPVJS0tT79695eXlpeDgYD355JO6dOmS3fGx5wkAAAAAAC6sKpy2M3jwYJ08eVIbN25UYWGh4uPjNWrUKK1cudLmuJEjR2r69OmWey8vL8s/FxUVqXfv3goNDdV3332nkydPKi4uTjVr1tTzzz9vV3wkTwAAAAAAQKXZv3+/1q9fr+3bt6t9+/aSpFdeeUW9evXS7NmzFRYWVupYLy8vhYaGlvhsw4YN2rdvnz7//HOFhIQoMjJSM2bM0MSJEzV16lTVqlWrzDFe++knAAAAAABQbmY3N6ddJpNJubm5VpfJZHIovpSUFPn5+VkSJ5IUExMjNzc3bd261ebYFStWKDAwUK1atVJiYqLy8vKs5r3pppsUEhJiaYuNjVVubq727t1rV4xlXnny8ssv2zWxJMXHx6tOnTp2jwMAAAAAAM5R7MSyneTkZE2bNs2qLSkpSVOnTi33nOnp6QoODrZqq1GjhgICApSenl7quPvvv1/h4eEKCwvTrl27NHHiRB08eFAffPCBZd7/TZxIstzbmrckZU6ePP7442rQoIHc3ct2PvQvv/yiPn36kDwBAAAAAMBFJCYmKiEhwarNaDSW2Pfpp5/WzJkzbc63f//+cscyatQoyz/fdNNNqlevnrp166YjR46oadOm5Z63JHbtefLDDz9clg0qDUkTAAAAAAAqn9mJp+QYjcZSkyV/NmHCBA0bNsxmnyZNmig0NFSZmZlW7ZcuXdKZM2dK3c+kJB07dpQkHT58WE2bNlVoaKi2bdtm1ScjI0OS7JpXsiN5kpSUJG9v7zJPPGnSJAUEBNgVDAAAAAAAcC6zm6FSPjcoKEhBQUFX7BcdHa3s7Gzt2LFDUVFRkqQvvvhCxcXFloRIWaSmpkqS6tWrZ5n3ueeeU2ZmpmUhyMaNG+Xj46OWLVva9S5lTj8lJSVZHflzJYmJifLz87MrGAAAAAAAUL20aNFCPXr00MiRI7Vt2zZ9++23Gjt2rAYNGmQ5aefXX3/VDTfcYFlJcuTIEc2YMUM7duzQsWPH9K9//UtxcXH6y1/+otatW0uSunfvrpYtW+qBBx7Qjz/+qM8++0yTJ0/WI488UubVM3/gqGJUaWazWT98+I4OfLlJprwLCr3uBt0WN0K+ofVKHVNcXKwdH76jQylfKy8nW7X9AnT9rV3Uru8AGQz/zcie/e2Etr6zQicP7lNxUbH86zfQnWMnqE7dwKvxagCqqT6DB+meEcPkHxSonw8c1BvTk/XTrj0l9m3UrKkeePwRXXdjS4U0qK83n5upNUv/adXHs7aX4h4fq+g7u8mvboCO7DugN//+D/20274d5gGgom3fe1SLPvpae478qlNnz+m1iUMU09G+vxkGULJi98pZeWKPFStWaOzYserWrZvc3Nw0YMAAq4NrCgsLdfDgQctpOrVq1dLnn3+uuXPn6sKFC2rYsKEGDBigyZMnW8a4u7tr7dq1euihhxQdHa3atWtr6NChmj59ut3xOTV5sn//fvXu3Vs///yzM6cFSvXjJx9pz8ZPdfvIR1QnKFjbP1itdS8+p3ufm6MapZzZnbpujfZt3qiuIx5RQP0GOnXsZ21Z9LpqeXnppjt7SZJyMtP10XNTdMNf7lD7u+9VTU9Pnf31hGrUrHk1Xw9ANfOXXrEaNelJvTJlhg7+uEv9hz6gvy9+UyO736WcM2cu6+/h6aH0X07om083aNSkp0qcc9xz0xRxfTPNfnKSsjIydUe/Pnp+2Vsa3bO/sjIySxwDAJUhz1Sg5hGhGnBHlMa+sKKywwFcSmWV7dgjICBAK1euLPV5RESEzGaz5b5hw4b68ssvrzhveHi4PvnkE4fjc96uMZIKCgp0/PhxZ04JlMpsNmv3hk/Uru9fFdHuZtVtGK7bR45V3tmzOrZze6njMg7/pPC27RUe2U51goLV5OZOanBja2X+fNjSZ/t7q9SodVt1GjhEgeGN5Rscqoi27eXp43s1Xg1ANXX3g3H6dPX72vj+GqUd/lmvTJkuU36+ut9zd4n9f9q9V4tmztGX69arsKDgsue1jEbdGhujRS/M0Z7tO3Qy7ReteOUN/Xb8F/W+f2BFvw4A2KVLu+Yaf3933dnpxsoOBQAuY9fKkz8fR/Rnp06dcigYwB7nTmUqLydb9Vu2trQZvbwU3LSZMo78pGadbilxXEiz67V/yyZlp/8mv9AwZaUdU/qhg4oeFCdJMhcXK23XTrXp2VfrZj+n08ePyicoWJG9+6txVIer8m4Aqp8aNWvouhtb6p35iyxtZrNZqd99rxZt25RrTvca7nKvUUOFJuvESsHFi7oxqq1D8QIAgKqjKqw8udbZlTyZN2+eIiMj5ePjU+Lz8+fPOyUooCzycrIlSZ6+1qtBPH18Lc9K0rZ3fxXm52t14ni5ubmpuLhYHQYM0nWdb5Mk5efmqvDiRaWu+0g3Dxiojn8brF92p2rDqy/qrolJCruB2lsAzufj7y/3GjV09nSWVfvZrCw1aNq4XHPmX8jTvp2puu+R0Uo78rOyT2epS59euqFtG508nuaMsAEAQBVgdq/sCKo+u5InzZo10/jx4zVkyJASn6emplqOFbLFZDLJZDJZtV0qKCh1jwpAkg5997W+WrbAct9zfGK55jmyLUWHvv9G3UY/Jv/6DZWVdkzfrVwqLz9/Nb+1q8zmYklSRLv2ah3bR5IUGB6hjMMHtW/zBpInAKqU2U8manzyDK349gsVXbqkw3v368u1n6pZK37LAAAAysqu5En79u21Y8eOUpMnBoPBagOX0iQnJ2vatGlWbd0fHK3YEQ/ZEw6qmfC27XVP0+ss90WXCiVJ+Tk5qu3nb2nPz81R3UYRpc7z/Tv/VGSvfpaynroNG+l81imlrl2j5rd2lUcdH7m5u8s/rIHVOL+w+kr/6aAT3wgA/iv37FkVXbok/8C6Vu3+devq7KmsUkZd2cm0E3pqcLyMnp7y8q6ts6dO6+m5s5T+ywlHQwYAAFUEZTuOs2vD2BdffFGPP/54qc/btGmj4uLiK86TmJionJwcq6tb3HB7QkE1VMvTU74hoZbLP6yBvHz99Ou+3ZY+Bfl5yjxyWCFNry91nksmkwxu1v/XN7i5WRJ/7jVqKKhxU2Wf/M2qT076SdUJ5JhiABXjUuElHdq7T5HRHS1tBoNBkZ07af+/f3R4flN+vs6eOi1vHx9F3dZZ33++2eE5AQBAFeHmxKuasmvlSWhoqFM+1Gg0ymg0WgdCyQ7sZDAYdFP3Xtr58QfyDa2nOoHB+uGDVfLy91dEu5st/T6eOV2NozqoVUwPSVJ4ZJT+/fEH8g4IVED9Bjqddky7Plur5rfdbhnTpmdfff76S6rXvIXCWrTSL7tTdTx1h+56eurVfk0A1ciHi5drwgvP6dCevTq4a7f6D3tARk9PbXx/jSRpwgvPKSsjU0tfnCfp901mGzVr+p9/rqm6IcFq0qK58i/k6WTaL5Kkdrd2lsFg0ImjxxQW3kjDJyboxM9HteE/cwLAteJCvklp6f9daXci84z2H/1Nvt5eCgvyq7zAAEB2Jk+Aa02bXv1UaDLpqyVvqiAvT6HX36BeEyZZJeNyMzN08Vyu5f6WIQ9q+wer9c3bC5Wfm6PafgFq0fVORfW7x9KncVQH3TZ0pP69bo2+XbFEfqFh6j52gupdf8NVfT8A1ctXn3wm34AADRn3iAKCAnVk/wE9O3yMsrN+/8NEcFg9q/LYgOBgvfav9yz394yI1z0j4rVr63ZNHPKgJKl2nTqKf2KcAkNDdC47R9989rmWzXlZRZcuXd2XA4Ar2HPkV8VNWWi5T17yiSTp7tvb6R+P3lPaMABlwYaxDjOYy7JJiaSAgAD99NNPCixj2UKjRo309ddfKzw8vEz956Q4viQZAK4lG+NK3h8KAKqqT9dMrewQAMC5bhxQ2RFcFZ2mO+/P299PaeO0uaqSMq88yc7O1qeffirfPx0LW5qsrCwVFRWVOzAAAAAAAIBrgV1lO0OHDq2oOAAAAAAAQEWoxhu9OkuZkydlOUUHAAAAAADA1ZR7w9hNmzZp06ZNyszMtEqsGAwGLVq0yCnBAQAAAAAAVLZyJU+mTZum6dOnq3379qpXr54MBoOz4wIAAAAAAE5goGzHYeVKnsyfP19Lly7VAw884Ox4AAAAAACAExncynTILmwoV/6poKBAnTt3dnYsAAAAAAAA15xyJU9GjBihlStXOjsWAAAAAADgZAY3513VVbnKdi5evKgFCxbo888/V+vWrVWzZk2r53PmzHFKcAAAAAAAwDFu7pUdQdVXruTJrl27FBkZKUnas2eP1TM2jwUAAAAAAK6kXMmTzZs3OzsOAAAAAABQAdyqcbmNs5QreQIAAAAAAKoGTttxHPknAAAAAAAAG1h5AgAAAACAC6Nsx3EkTwAAAAAAcGEkTxzHVwgAAAAAAGADK08AAAAAAHBhrDxxHMkTAAAAAABcGMkTx/EVAgAAAAAA2MDKEwAAAAAAXBgrTxxH8gQAAAAAABfm7mau7BCqPPJPAAAAAAAANrDyBAAAAAAAF0bZjuNIngAAAAAA4MJInjiOrxAAAAAAAMAGkicAAAAAAAA2ULYDAAAAAIALc2fZhMP4CgEAAAAAAGxg5QkAAAAAAC7MzVDZEVR9JE8AAAAAAHBhlO04jq8QAAAAAADABpInAAAAAAC4MDc3510V5cyZMxo8eLB8fHzk5+en4cOH6/z586X2P3bsmAwGQ4nXu+++a+lX0vNVq1bZHR9lOwAAAAAAuLCqULYzePBgnTx5Uhs3blRhYaHi4+M1atQorVy5ssT+DRs21MmTJ63aFixYoFmzZqlnz55W7UuWLFGPHj0s935+fnbHR/IEAAAAAABUmv3792v9+vXavn272rdvL0l65ZVX1KtXL82ePVthYWGXjXF3d1doaKhV24cffqh7771X3t7eVu1+fn6X9bVXFcg/AQAAAACA8nJ3c95VEVJSUuTn52dJnEhSTEyM3NzctHXr1jLNsWPHDqWmpmr48OGXPXvkkUcUGBioDh06aPHixTKbzXbHyMoTAAAAAABcmDOTHiaTSSaTyarNaDTKaDSWe8709HQFBwdbtdWoUUMBAQFKT08v0xyLFi1SixYt1LlzZ6v26dOn64477pCXl5c2bNighx9+WOfPn9djjz1mV4ysPAEAAAAAAGWSnJwsX19fqys5ObnEvk8//XSpm7r+cR04cMDhmPLz87Vy5coSV508++yzuuWWW9S2bVtNnDhRTz31lGbNmmX3Z7DyBAAAAAAAF+bMU3ISExOVkJBg1VbaqpMJEyZo2LBhNudr0qSJQkNDlZmZadV+6dIlnTlzpkx7lbz33nvKy8tTXFzcFft27NhRM2bMkMlksmu1DMkTAAAAAABcmLvBeXPZU6ITFBSkoKCgK/aLjo5Wdna2duzYoaioKEnSF198oeLiYnXs2PGK4xctWqS+ffuW6bNSU1Pl7+9vd5kRyRMAAAAAAFBpWrRooR49emjkyJGaP3++CgsLNXbsWA0aNMhy0s6vv/6qbt26afny5erQoYNl7OHDh/XVV1/pk08+uWzejz/+WBkZGerUqZM8PDy0ceNGPf/883riiSfsjpHkCQAAAAAALqyiTslxphUrVmjs2LHq1q2b3NzcNGDAAL388suW54WFhTp48KDy8vKsxi1evFgNGjRQ9+7dL5uzZs2aeu211zR+/HiZzWY1a9ZMc+bM0ciRI+2Oz2Auzxk9FWBOyo+VHQIAONXGuCGVHQIAONWna6ZWdggA4Fw3DqjsCK6Kcet2Om2ueb3bOW2uqqQK5J8AAAAAAAAqD2U7AAAAAAC4sBpuTtwxtpoieQIAAAAAgAurCnueXOv4CgEAAAAAAGxg5QkAAAAAAC7Mnaodh5E8AQAAAADAhVG24zi+QgAAAAAAABtIngAAAAAAANhA2Q4AAAAAAC6Msh3H8RUCAAAAAADYwMoTAAAAAABcmLsbx+04iuQJAAAAAAAujLIdx/EVAgAAAAAA2MDKEwAAAAAAXJg7VTsOI3kCAAAAAIALY88Tx1G2AwAAAAAAYAMrTwAAAAAAcGFsGOs4g9lsNld2EMDVYjKZlJycrMTERBmNxsoOBwAcxu8aAFfD7xqAaxHJE1Qrubm58vX1VU5Ojnx8fCo7HABwGL9rAFwNv2sArkUs3gEAAAAAALCB5AkAAAAAAIANJE8AAAAAAABsIHmCasVoNCopKYnNxwC4DH7XALgaftcAXIvYMBYAAAAAAMAGVp4AAAAAAADYQPIEAAAAAADABpInAAAAAAAANpA8gUuLiIiQwWCQwWBQdnZ2mcctXbrUMu7xxx+vsPgAwF7l/V2bOnWqZdzcuXMrLD4AsNcfv01+fn52jeN3DcDVRPIELm/69Ok6efKkfH19JUkXL17UsGHDdNNNN6lGjRrq37//ZWMGDhyokydPKjo6+ipHCwBX9ufftS1btqhfv36qV6+eateurcjISK1YscJqzBNPPKGTJ0+qQYMGlREyANi0ZMkS/fTTT5b7kydP6v7779f1118vNze3Ev8yi981AFcTyRO4vDp16ig0NFQGg0GSVFRUJE9PTz322GOKiYkpcYynp6dCQ0NVq1atqxkqAJTJn3/XvvvuO7Vu3Vrvv/++du3apfj4eMXFxWnt2rWWMd7e3goNDZW7u3tlhQ0ApfLz81NwcLDl3mQyKSgoSJMnT1abNm1KHMPvGoCrqUZlBwA4omvXrmrVqpUk6e2331bNmjX10EMPafr06ZY/VPxZ7dq19cYbb0iSvv32W7uWvQNARSvP79qkSZOs7seNG6cNGzbogw8+UJ8+fSo8ZgCwpWvXrmrdurU8PDy0cOFC1apVS2PGjNHUqVNLHRMREaF58+ZJkhYvXnyVIgWA0rHyBFXesmXLVKNGDW3btk3z5s3TnDlztHDhwsoOCwDKzRm/azk5OQoICKigCAHAPsuWLVPt2rW1detWvfDCC5o+fbo2btxY2WEBQJmx8gRVXsOGDfXSSy/JYDCoefPm2r17t1566SWNHDmyskMDgHJx9HftnXfe0fbt2/Xmm29WcKQAUDatW7dWUlKSJOm6667Tq6++qk2bNunOO++s5MgAoGxYeYIqr1OnTlZL2aOjo3Xo0CEVFRVVYlQAUH6O/K5t3rxZ8fHxeuutt3TjjTdWZJgAUGatW7e2uq9Xr54yMzMrKRoAsB/JEwAAXMSXX36pu+66Sy+99JLi4uIqOxwAsKhZs6bVvcFgUHFxcSVFAwD2I3mCKm/r1q1W999//72uu+46dl4HUGWV53dty5Yt6t27t2bOnKlRo0ZVdIgAAADVCnueoMpLS0tTQkKCRo8erZ07d+qVV17Riy++aHPMvn37VFBQoDNnzujcuXNKTU2VJEVGRlZ8wABwBfb+rm3evFl9+vTRuHHjNGDAAKWnp0uSatWqxaaxAKqsP/777Pz58zp16pRSU1NVq1YttWzZsnIDA1AtkTxBlRcXF6f8/Hx16NBB7u7uGjdu3BX/1rVXr146fvy45b5t27aSJLPZXKGxAkBZ2Pu7tmzZMuXl5Sk5OVnJycmW9i5dumjLli1XIWIAcL4//vtMknbs2KGVK1cqPDxcx44dq7ygAFRbJE9Q5dWsWVNz587VG2+8UeYx/EsXwLXM3t+1pUuXaunSpRUbFACUU0lJ3DVr1lxxHH+pBeBawp4ncHkTJ06Ut7e3cnJyyjxmxYoV8vb21tdff12BkQFA+ZTnd+3555+Xt7e30tLSKjAyACif++67Tw0aNLBrDL9rAK4mVp7ApX355ZcqLCyUJNWpU6fM4/r27auOHTtKkvz8/CoiNAAol/L+ro0ZM0b33nuvJCkoKKhCYgOA8jh06JAk2b3ZP79rAK4mg5n1cAAAAAAAAKWibAcAAAAAAMAGkicAAAAAAAA2kDwBAAAAAACwgeQJAAAAAACADSRPAAAAAAAAbCB5AgAAAAAAYAPJEwAAAAAAABtIngAAAAAAANhA8gQAAAAAAMCG/weAyEmeOcThwgAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 1500x500 with 2 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(15, 5))\n",
        "sns.heatmap(df[['p[1]','p[2]','n[1]']].corr(),annot=True, center=0,ax=axes)\n",
        "\n",
        "axes.set_title('Correlations')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a4aa96d4",
      "metadata": {
        "id": "a4aa96d4"
      },
      "source": [
        "In this problem, we've assumed that the amount of space we have available for the products is 200 units. In retail, this could be the amount of warehouse space, or for ticketing this could represent the number of seats available."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6e9c6157",
      "metadata": {
        "id": "6e9c6157"
      },
      "source": [
        "### Building regressors to predict sales\n",
        "\n",
        "The prices for each category item will be used to predict the number of Category 1 items sold. Here we build a regression model to form this relationship which will later be used as part of the optimization model."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "a70d5f2c",
      "metadata": {
        "id": "a70d5f2c"
      },
      "outputs": [],
      "source": [
        "from sklearn.compose import make_column_transformer\n",
        "from sklearn.linear_model import LinearRegression\n",
        "from sklearn.pipeline import make_pipeline\n",
        "from sklearn.metrics import r2_score\n",
        "from sklearn.model_selection import train_test_split      #importing scikit-learn's function for data splitting\n",
        "from sklearn.ensemble import GradientBoostingRegressor    #importing scikit-learn's gradient booster regressor function\n",
        "from sklearn.model_selection import cross_validate        #improting scikit-learn's cross validation function"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "3095da93",
      "metadata": {
        "id": "3095da93"
      },
      "outputs": [],
      "source": [
        "X = df[[\"p[1]\",\"p[2]\"]]\n",
        "y = df[\"n[1]\"]\n",
        "# Split the data for training and testing\n",
        "X_train, X_test, y_train, y_test = train_test_split(\n",
        "    X, y, train_size=0.75, random_state=1\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f8d713c2",
      "metadata": {
        "id": "f8d713c2"
      },
      "source": [
        "First we'll start with a linear regression model."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "33f4eaec",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "33f4eaec",
        "outputId": "d02c9203-4bef-4a4e-cb16-bb948a8fca00"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(array([0.80665368, 0.81004987, 0.7941077 , 0.79780172, 0.79495142]),\n",
              " array([0.77650521, 0.74227784, 0.82105606, 0.81102818, 0.82166193]))"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "linear_regressor = make_pipeline(LinearRegression())\n",
        "linear_regressor.fit(X_train, y_train)\n",
        "linear_regression_validation = cross_validate(linear_regressor, X_train, y_train, cv=5, return_train_score=True, return_estimator=True)\n",
        "\n",
        "linear_regression_validation['train_score'],linear_regression_validation['test_score']"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cc9dd59b",
      "metadata": {
        "id": "cc9dd59b"
      },
      "source": [
        "Let's try a gradient boosting model as well."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "9b8eb814",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "9b8eb814",
        "outputId": "ca16f36f-2399-4e29-c3b3-9646cdb22a38"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(array([0.71068886, 0.71262322, 0.70264476, 0.70564959, 0.7008412 ]),\n",
              " array([0.65706085, 0.69348053, 0.66690978, 0.67664668, 0.67992755]))"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "from sklearn.ensemble import GradientBoostingRegressor\n",
        "xgb_regressor = make_pipeline(GradientBoostingRegressor(n_estimators=10))\n",
        "xgb_regressor.fit(X_train, y_train)\n",
        "xgb_regressor_validation = cross_validate(xgb_regressor, X_train, y_train, cv=5, return_train_score=True, return_estimator=True)\n",
        "\n",
        "xgb_regressor_validation['train_score'], xgb_regressor_validation['test_score']\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "mTquaGiJF2pO",
      "metadata": {
        "id": "mTquaGiJF2pO"
      },
      "source": [
        "## Price optimization model with competing products\n",
        "\n",
        "Our problem is to:\n",
        "1.\tDetermine the number of each category of product to make available given the overall restriction of what we can offer.\n",
        "2.\tWe are also instructed to make sure there are a minimum number of each category made available as well as a minimum and maximum price for each category.\n",
        "3.\tLastly, the product categories should be decreasing in price, meaning Category 1 should be the most expensive, and so on. Specifically, we must make sure there is at least a $50 gap between categories, but no more than $100.\n",
        "\n",
        "With the predictive part in place, it's time to build the optimization model. The model is formulated (i.e. the mathematical representation) for an unspecified number of categories, but the code will reflect that we have two categories of products in this problem. We start by setting some parameter values (not to be confused with ML hyperparameters) and initialize the optimization model.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "95b966b2",
      "metadata": {
        "id": "95b966b2"
      },
      "source": [
        "### Initialize model and set input parameters\n",
        "- $C$: Number of product categories\n",
        "- $N$: Total amount of space available\n",
        "- $\\lambda$: Price control parameter\n",
        "\n",
        "Here is the first mention of a price control parameter. It is fairly common in optimization modeling to add penalty terms to try and prevent undesirable outcomes. This is akin to using penalty terms in machine learning and applied statistics to prevent overfitting, with [Lasso](https://en.wikipedia.org/wiki/Lasso_(statistics)) and [Ridge](https://en.wikipedia.org/wiki/Ridge_regression) regression as a couple of common examples."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "67db9804",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "67db9804",
        "outputId": "2a1edc1c-0b23-4649-88da-53d34a5f2b03"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Restricted license - for non-production use only - expires 2025-11-24\n"
          ]
        }
      ],
      "source": [
        "#### Initialize the model\n",
        "m = gp.Model(\"price optimization\")\n",
        "\n",
        "products = [1,2]            #### Category 1 and Category 2\n",
        "N = 200                     #### limit on available space\n",
        "l = 0                       #### price control, we'll start this at 0"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ac5a875b",
      "metadata": {
        "id": "ac5a875b"
      },
      "source": [
        "### Decision variables\n",
        "- $p_c$: price per item in category $c = 1,2,\\dots, C$\n",
        "- $n_c$: number of items allocated to category, predicted using features $p_c$, $c = 1,2,\\dots, C$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "bbe94a3c",
      "metadata": {
        "id": "bbe94a3c"
      },
      "outputs": [],
      "source": [
        "p = m.addVars(products, name=\"p\")            #### price decision variables\n",
        "n = m.addVars(products, name=\"n\")            #### decision variable for number of items in each category"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "721614e0",
      "metadata": {
        "id": "721614e0"
      },
      "source": [
        "### Constraints"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8403f941",
      "metadata": {
        "id": "8403f941"
      },
      "source": [
        "We need to have a minimum number of each category available.\n",
        "\\begin{align*}\n",
        "n_c \\ge l_{n_c}\n",
        "\\end{align*}\n",
        "\n",
        "We also set lower and upper bounds on the prices.\n",
        "\\begin{align*}\n",
        "l_{n_c} \\le p_c \\le u_{p_c}\n",
        "\\end{align*}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "edec4432",
      "metadata": {
        "id": "edec4432"
      },
      "outputs": [],
      "source": [
        "min_items = {1:50,2:50}\n",
        "price_bounds = {1:[300,400], 2:[100,300]}\n",
        "m.addConstrs(n[c] >= min_items[c] for c in products)        #### we could hardcode 50 instead of min_items, but this is more flexible\n",
        "m.addConstr(p[1] == [300,400])                              #### this is a shorthand way to code 300 <= p[1] <= 400\n",
        "m.addConstr(p[2] == [100,300]);"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3bef9e0d",
      "metadata": {
        "id": "3bef9e0d"
      },
      "source": [
        "Another note: each of the above constraints can be addressed when defining the decision variables. Here is an example for the decision variable $n$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "e2065d83",
      "metadata": {
        "id": "e2065d83"
      },
      "outputs": [],
      "source": [
        "#price_lb = {1:300, 2:100}\n",
        "#price_ub = {1:400, 2:300}\n",
        "#p = m.addVars(products, lb = price_lb, ub = price_ub, name=\"p\")     #### each price is now bounded\n",
        "#n = m.addVars(products, lb = min_items, name=\"n\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0578e271",
      "metadata": {
        "id": "0578e271"
      },
      "source": [
        "In general, the number of items allocated must equal the total available space.\n",
        "\\begin{equation*}\n",
        "n_1 + n_2 + \\dots + n_C = \\sum_{c}n_c = N \\\\\n",
        "\\end{equation*}\n",
        "Note that this, along with the constraint on the minimum number available means we don't have to specify an upper bound for each $n_c$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "140e739f",
      "metadata": {
        "id": "140e739f"
      },
      "outputs": [],
      "source": [
        "m.addConstr(n.sum() ==  N);     #### remember we set N = 200 earlier"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b17a1e78",
      "metadata": {
        "id": "b17a1e78"
      },
      "source": [
        "The last set of constraints are for price ordering. This requires the subsequent category to be cost between $50 and $100 less than the previous.\n",
        "\\begin{equation*}\n",
        "50 \\le p_c - p_{c+1} \\le 100\n",
        "\\end{equation*}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "1428be22",
      "metadata": {
        "id": "1428be22"
      },
      "outputs": [],
      "source": [
        "m.addConstr(p[1]-p[2] == [50,100]);"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "51ee2204",
      "metadata": {
        "id": "51ee2204"
      },
      "source": [
        "### Objective function"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6eb69e54",
      "metadata": {
        "id": "6eb69e54"
      },
      "source": [
        "We want to maximize total revenue with the portion of total revenue coming from category $c$ being $p_cn_c$. This makes the total revenue $\\sum_{c} p_c n_c$. That is the first part of the objective. Earlier a price control parameter was introduced which is the second part of the objective. The lambda parameter captures the trade-off between the revenue and price-control pieces. This term penalizes the model from setting too high of prices since doing so could lose sales. Our model assumed we'll sell all of the items so having this penalty term can make this assumption more realistic. For reference, [here is a good source](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4565407).\n",
        "\n",
        "This term will be defined as $λ (\\sum_{c} p_c^2)$ for this problem. So, the complete objective is:\n",
        "\\begin{equation*}\n",
        "\\textrm{maximize} \\sum_{c} p_c n_c - λ (\\sum_{c} p_c^2)\n",
        "\\end{equation*}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "a2cfcabc",
      "metadata": {
        "id": "a2cfcabc"
      },
      "outputs": [],
      "source": [
        "revenue = gp.quicksum(p[c]*n[c] for c in products)          #### you could also use the more simple p.prod(n)\n",
        "penalty = l*(p[1]**2+p[2]**2)                               #### we used l as the lambda parameter earlier\n",
        "m.setObjective(revenue - penalty, sense = GRB.MAXIMIZE)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c707af60",
      "metadata": {
        "id": "c707af60"
      },
      "source": [
        "### Integrate the ML model\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "38d41039",
      "metadata": {
        "id": "38d41039"
      },
      "source": [
        "Right now, if we were to run the optimization, the solution would be to set the price for Category 1 to $400, Category 2 to $300, and sell 150 and 50 of each item, respectively. That's because we have yet to add in the relationship between price and demand that was derived from the ML model. To integrate the machine learning model into the optimization model, we'll use the Gurobi Machine Learning package. The magic happens using `add_predictor_constr` function."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "b113eda2",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "b113eda2",
        "outputId": "584d9f0e-ea69-41c9-e481-a45f0bc36dd6"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: gurobi-machinelearning in /usr/local/lib/python3.10/dist-packages (1.4.0)\n",
            "Requirement already satisfied: numpy>=1.22.0 in /usr/local/lib/python3.10/dist-packages (from gurobi-machinelearning) (1.25.2)\n",
            "Requirement already satisfied: gurobipy>=10.0.0 in /usr/local/lib/python3.10/dist-packages (from gurobi-machinelearning) (11.0.0)\n",
            "Requirement already satisfied: scipy>=1.9.3 in /usr/local/lib/python3.10/dist-packages (from gurobi-machinelearning) (1.11.4)\n"
          ]
        }
      ],
      "source": [
        "#### install the package and load the required function\n",
        "%pip install gurobi-machinelearning\n",
        "from gurobi_ml import add_predictor_constr"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "119a8f7e",
      "metadata": {
        "id": "119a8f7e"
      },
      "source": [
        "This additional package is useful when we have **decision variables** that are also **features** of a machine learning model. First, we need a data frame that contains these decision variables. It is important to make sure the indices of the data frame have the **same name** as the training data for the machine learning model."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "0b1ac8cc",
      "metadata": {
        "id": "0b1ac8cc"
      },
      "outputs": [],
      "source": [
        "m_feats = pd.DataFrame({\"p[1]\":[p[1]],\"p[2]\":[p[2]]})"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3263a396",
      "metadata": {
        "id": "3263a396"
      },
      "source": [
        "Adding the predictive model to the optimization model requires specifying the model we want to use `(m)`, regression object `(xgb_regressor)`, feature data frame `(m_feats)`, and the output decision variable `(n[1])`. Remember `n[2]` is **NOT** the output of the regression. We can then print the number of variables and constraints added to the model using `print_stats`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "1fe6c59f",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1fe6c59f",
        "outputId": "5aa94f80-e974-4e5b-de19-7c82525e5f75"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Model for pipe:\n",
            "90 variables\n",
            "11 constraints\n",
            "246 general constraints\n",
            "Input has shape (1, 2)\n",
            "Output has shape (1, 1)\n",
            "\n",
            "Pipeline has 1 steps:\n",
            "\n",
            "--------------------------------------------------------------------------------\n",
            "Step            Output Shape    Variables              Constraints              \n",
            "                                                Linear    Quadratic      General\n",
            "================================================================================\n",
            "gbtree_reg            (1, 1)           90           11            0          246\n",
            "\n",
            "--------------------------------------------------------------------------------\n"
          ]
        }
      ],
      "source": [
        "pred_constr = add_predictor_constr(m, xgb_regressor, m_feats, n[1])\n",
        "pred_constr.print_stats()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "53ce50f6",
      "metadata": {
        "id": "53ce50f6"
      },
      "source": [
        "### Solve the optimization and get the solution"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ec80341d",
      "metadata": {
        "id": "ec80341d"
      },
      "source": [
        "Since this is a quadratic, non-convex problem we set the `NonConvex` parameter to 2. See the [documentation](https://www.gurobi.com/documentation/current/refman/nonconvex.html) for more information. We'll also print out the optimal solution."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "d04d5e5b",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "d04d5e5b",
        "outputId": "7cbae0da-917f-45ec-9fe1-3af04c942d08"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Set parameter NonConvex to value 2\n",
            "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.3 LTS\")\n",
            "\n",
            "CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n",
            "Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n",
            "\n",
            "Optimize a model with 17 rows, 97 columns and 102 nonzeros\n",
            "Model fingerprint: 0x2255329b\n",
            "Model has 2 quadratic objective terms\n",
            "Model has 246 general constraints\n",
            "Variable types: 17 continuous, 80 integer (80 binary)\n",
            "Coefficient statistics:\n",
            "  Matrix range     [1e-01, 1e+00]\n",
            "  Objective range  [0e+00, 0e+00]\n",
            "  QObjective range [2e+00, 2e+00]\n",
            "  Bounds range     [1e+00, 2e+02]\n",
            "  RHS range        [1e+00, 4e+02]\n",
            "  GenCon rhs range [3e-01, 4e+02]\n",
            "  GenCon coe range [1e+00, 1e+00]\n",
            "Presolve added 44 rows and 0 columns\n",
            "Presolve removed 0 rows and 15 columns\n",
            "Presolve time: 0.10s\n",
            "Presolved: 66 rows, 85 columns, 383 nonzeros\n",
            "Presolved model has 2 bilinear constraint(s)\n",
            "\n",
            "Solving non-convex MIQCP\n",
            "\n",
            "Variable types: 18 continuous, 67 integer (67 binary)\n",
            "Found heuristic solution: objective 57353.074762\n",
            "\n",
            "Root relaxation: objective 6.842235e+04, 81 iterations, 0.00 seconds (0.00 work units)\n",
            "\n",
            "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
            " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
            "\n",
            "     0     0 68422.3454    0   21 57353.0748 68422.3454  19.3%     -    0s\n",
            "     0     0 68380.4264    0   21 57353.0748 68380.4264  19.2%     -    0s\n",
            "H    0     0                    65051.531333 68380.4264  5.12%     -    0s\n",
            "H    0     0                    65717.039440 68380.4264  4.05%     -    0s\n",
            "     0     0 67521.7136    0   30 65717.0394 67521.7136  2.75%     -    0s\n",
            "     0     0 67440.6730    0   34 65717.0394 67440.6730  2.62%     -    0s\n",
            "     0     0 67415.9059    0   33 65717.0394 67415.9059  2.59%     -    0s\n",
            "     0     0 67407.6626    0   32 65717.0394 67407.6626  2.57%     -    0s\n",
            "     0     0 67397.1289    0   32 65717.0394 67397.1289  2.56%     -    0s\n",
            "     0     0 67367.1673    0   26 65717.0394 67367.1673  2.51%     -    0s\n",
            "H    0     0                    65919.660132 67367.1673  2.20%     -    0s\n",
            "     0     0 67325.6782    0   24 65919.6601 67325.6782  2.13%     -    0s\n",
            "     0     0 67027.1365    0   27 65919.6601 67027.1365  1.68%     -    0s\n",
            "     0     0 66962.1686    0   21 65919.6601 66962.1686  1.58%     -    0s\n",
            "     0     0 66754.1535    0   21 65919.6601 66754.1535  1.27%     -    0s\n",
            "     0     0 66724.4366    0   21 65919.6601 66724.4366  1.22%     -    0s\n",
            "     0     0 66357.6817    0   22 65919.6601 66357.6817  0.66%     -    0s\n",
            "     0     0 66314.3792    0   17 65919.6601 66314.3792  0.60%     -    0s\n",
            "     0     0 66212.7588    0   15 65919.6601 66212.7588  0.44%     -    0s\n",
            "     0     0 66096.5119    0   21 65919.6601 66096.5119  0.27%     -    0s\n",
            "     0     0 66096.5119    0   22 65919.6601 66096.5119  0.27%     -    0s\n",
            "     0     0 66096.5119    0   21 65919.6601 66096.5119  0.27%     -    0s\n",
            "     0     0 66096.5119    0   23 65919.6601 66096.5119  0.27%     -    0s\n",
            "     0     0 66092.0244    0   21 65919.6601 66092.0244  0.26%     -    0s\n",
            "     0     0 66092.0244    0   23 65919.6601 66092.0244  0.26%     -    0s\n",
            "     0     0 66092.0244    0   22 65919.6601 66092.0244  0.26%     -    0s\n",
            "     0     0 66092.0244    0   21 65919.6601 66092.0244  0.26%     -    0s\n",
            "     0     0 66092.0244    0   21 65919.6601 66092.0244  0.26%     -    0s\n",
            "     0     0 66092.0244    0   23 65919.6601 66092.0244  0.26%     -    0s\n",
            "     0     0 66092.0244    0   19 65919.6601 66092.0244  0.26%     -    0s\n",
            "     0     0 66087.8039    0   21 65919.6601 66087.8039  0.26%     -    0s\n",
            "     0     0 66078.8922    0   17 65919.6601 66078.8922  0.24%     -    0s\n",
            "     0     0 66074.6109    0   20 65919.6601 66074.6109  0.24%     -    0s\n",
            "     0     0 65978.7306    0   15 65919.6601 65978.7306  0.09%     -    0s\n",
            "     0     0 65955.6109    0   23 65919.6601 65955.6109  0.05%     -    0s\n",
            "     0     0 65941.1212    0   23 65919.6601 65941.1212  0.03%     -    0s\n",
            "     0     0 65933.9387    0   23 65919.6601 65933.9387  0.02%     -    0s\n",
            "     0     0 65933.9387    0   22 65919.6601 65933.9387  0.02%     -    0s\n",
            "     0     0 65931.8301    0    1 65919.6601 65931.8301  0.02%     -    0s\n",
            "H    0     0                    65931.830113 65931.8301  0.00%     -    0s\n",
            "     0     0 65931.8301    0    1 65931.8301 65931.8301  0.00%     -    0s\n",
            "\n",
            "Cutting planes:\n",
            "  Cover: 3\n",
            "  Implied bound: 2\n",
            "  Clique: 6\n",
            "  MIR: 2\n",
            "  GUB cover: 1\n",
            "  Relax-and-lift: 4\n",
            "\n",
            "Explored 1 nodes (377 simplex iterations) in 0.62 seconds (0.05 work units)\n",
            "Thread count was 2 (of 2 available processors)\n",
            "\n",
            "Solution count 5: 65931.8 65919.7 65717 ... 57353.1\n",
            "\n",
            "Optimal solution found (tolerance 1.00e-04)\n",
            "Best objective 6.593183011274e+04, best bound 6.593183011274e+04, gap 0.0000%\n"
          ]
        }
      ],
      "source": [
        "m.Params.NonConvex = 2\n",
        "m.optimize()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "ec56f1da",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ec56f1da",
        "outputId": "231a9f5b-2d76-48b2-899b-8f9db2232b94"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "Optimal price for the two categories:\n",
            " 364.97 300.0\n",
            "\n",
            "Optimal number of space assigned to the two categories:\n",
            " 91 109\n",
            "\n",
            "Total revenue:\n",
            " 65931.83\n"
          ]
        }
      ],
      "source": [
        "print(\"\\nOptimal price for the two categories:\\n\",round(p[1].X,2),round(p[2].X,2))\n",
        "print(\"\\nOptimal number of space assigned to the two categories:\\n\",round(n[1].X), round(n[2].X))\n",
        "print(\"\\nTotal revenue:\\n\",round(revenue.getValue(),2))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5e9e2bbd",
      "metadata": {
        "id": "5e9e2bbd"
      },
      "source": [
        "## Follow-up questions\n",
        "At this point the λ value is set to zero, meaning we are not penalizing the price. One way to see the sensitivity of the prices and number of items set to each category is solving the model for multiple values of λ. This can be done fairly easy in `gurobipy`. Check out the documentation for how Gurobi can handle [multiple scenarios]('https://www.gurobi.com/documentation/current/refman/multiple_scenarios.html').\n",
        "\n",
        "We also only used the trained `xgb_regressor` for the optimization to show the number of variables and constraints added to the model. Test the `linear_regressor` that was trained as well. You can do that by re-running the cells or by adding to the code to [remove the previous predictor constraints](https://gurobi-machinelearning.readthedocs.io/en/stable/auto_generated/gurobi_ml.sklearn.linear_regression.LinearRegressionConstr.html#gurobi_ml.sklearn.linear_regression.LinearRegressionConstr.remove) and adding new ones."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e0888fb6",
      "metadata": {
        "id": "e0888fb6"
      },
      "source": [
        "Copyright © 2024 Gurobi Optimization, LLC"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}
